home *** CD-ROM | disk | FTP | other *** search
/ Libris Britannia 4 / science library(b).zip / science library(b) / MATH / SPECTR20.ZIP / COLOR2.FOR < prev    next >
Text File  |  1992-04-23  |  9KB  |  342 lines
  1. *    COLOR2.FOR
  2.  
  3. *    Program to produce 2 channels of colored noise for spectrum test.
  4.  
  5. *    David E. Hess
  6. *    Fluid Flow Group - Process Measurements Division
  7. *    Chemical Science and Technology Laboratory
  8. *    National Institute of Standards and Technology
  9. *    April 15, 1992
  10.  
  11.     IMPLICIT    REAL*4 (A-H,O-Z), INTEGER*2 (I-N)
  12.     PARAMETER     (NUMO=2,NSMAX=20,NMAX=8192)
  13.     INTEGER*2     NDATA[ALLOCATABLE,HUGE](:)
  14.     INTEGER*2    GAIN(0:7)
  15.     INTEGER*4     IRSIZE,N,J,ITST,IDELTMS
  16.     REAL*4         X(NSMAX+1,5),X2(NSMAX+1,5)
  17.     REAL*4         A(NSMAX),B(NSMAX),C(NSMAX)
  18.     REAL*4         D(NSMAX),E(NSMAX),GRAF(2,20)
  19.     REAL*4         A2(NSMAX),B2(NSMAX),C2(NSMAX)
  20.     REAL*4         D2(NSMAX),E2(NSMAX),GRAF2(2,20)
  21.     CHARACTER*1    FIRST
  22.     CHARACTER*4     FLNM
  23.     CHARACTER*8     FILNAM
  24.  
  25. *    This routine uses a random number generator to produce uniform
  26. *    deviates between 0 and 1.  These are then transformed to a 
  27. *    sequence having a power density = 1 with power concentrated
  28. *    between f1 and f2 Hz for a sampling interval of dt secs. This 
  29. *    process is then repeated to yield a second channel.
  30.  
  31. *    Initialization.
  32.  
  33.     ICHANS=2
  34.     GAIN=0
  35.     VOFST=20.0/4096.0
  36.  
  37.     WRITE (*,'(/1X,A)') 'Channel 1'
  38.     WRITE (*,'(1X,A)')  '---------'
  39.     WRITE (*,'(/1X,A\)') 'Enter the low frequency cutoff : '
  40.     READ (*,*) F1
  41.     WRITE (*,'(/1X,A\)') 'Enter the high frequency cutoff : '
  42.     READ (*,*) F2
  43.     WRITE (*,'(/1X,A)') 'Channel 2'
  44.     WRITE (*,'(1X,A)')  '---------'
  45.     WRITE (*,'(/1X,A\)') 'Enter the low frequency cutoff : '
  46.     READ (*,*) F12
  47.     WRITE (*,'(/1X,A\)') 'Enter the high frequency cutoff : '
  48.     READ (*,*) F22
  49.  
  50. *    Get the sampling interval.
  51.  
  52. 10    WRITE (*,'(/1X,A\)') 'Enter delta T secs. (1.0/Sample Rate) : '
  53.     READ (*,*) DELT
  54.       XTST=SQRT(6.0/DELT)*0.5
  55.       ITST=INT4(XTST/VOFST)
  56.       IF (ITST .GT. 32767) THEN
  57.         WRITE (*,'(//1X,A,F7.6,A)')
  58.      +        'This choice of delta t = ',DELT,' ,'
  59.         WRITE (*,'(1X,A,F6.2,A)')
  60.      +        'leads to a maximum data value of ',XTST,' ,'
  61.         WRITE (*,'(1X,A,A,I6,A)')
  62.      +        'which, when converted to an integer,',
  63.      +        ' has a value of ',ITST,' .'
  64.         WRITE (*,'(1X,A/)')
  65.      +        'This cannot be stored in the Integer*2 array NDATA.'
  66.         GO TO 10
  67.       ENDIF
  68.  
  69. *    Get order of filter.
  70.  
  71. 20    WRITE (*,'(/1X,A\)') 'Enter # of filter sections to cascade : '
  72.     READ (*,*) NS
  73.     IF (NS .GT. NSMAX) THEN
  74.       WRITE (*,'(1X,A,I2)') 'Reduce the # of sections to ',NSMAX
  75.       GO TO 20
  76.     ENDIF
  77.  
  78. *    Get number of points per record.
  79.  
  80. 30    WRITE (*,'(/1X,A\)')
  81.      +      'Enter # of points per channel per record (N) : '
  82.     READ (*,*) N
  83.     IF (N .GT. NMAX) THEN
  84.       WRITE (*,'(1X,A,I5)') 'Reduce the # of points to ',NMAX
  85.       GO TO 30
  86.     ENDIF
  87.  
  88. *    Allocate space for NDATA array.
  89.  
  90.     ALLOCATE (NDATA(2*N), STAT=IERR)
  91.     IF (IERR .NE. 0)
  92.      +          STOP 'Not enough storage for data.  Aborting ...'
  93.  
  94. *    Get number of records in file.
  95.  
  96.     WRITE (*,'(/1X,A\)') 'Enter # of records (NUMREC) : '
  97.     READ (*,*) NUMREC
  98.  
  99. *    Get seed value for random number generator.
  100.  
  101.     WRITE (*,'(/1X,A\)') 'Enter a negative integer : '
  102.     READ (*,*) IDUM
  103.  
  104. *    Get output file name.
  105.  
  106. 40    WRITE (*,'(/1X,A\)') 'Enter output file name (4 chars) : '
  107.     READ (*,'(A)') FLNM
  108.     WRITE (*,'( )')
  109.  
  110. *    Convert to uppercase and check first character alphabetic.
  111.  
  112.     DO J=4,1,-1
  113.       FIRST=FLNM(J:J)
  114.       IF (ICHAR(FIRST) .GE. 97 .AND. ICHAR(FIRST) .LE. 122) THEN
  115.         IHOLD=ICHAR(FIRST)-32
  116.         FIRST=CHAR(IHOLD)
  117.         FLNM(J:J)=FIRST
  118.       ENDIF
  119.     ENDDO
  120.     IF (ICHAR(FIRST) .LT. 65 .OR. ICHAR(FIRST) .GT. 90) THEN
  121.       WRITE (*,'(/1X,A,A,A/1X,A,A,A/1X,A)') 
  122.      +      'Filename ',FLNM,' began with',
  123.      +      'the nonalphabetic character ',FIRST,'.',
  124.      +      'Re-enter the filename correctly.'
  125.       GO TO 40
  126.     ENDIF
  127.  
  128.     FILNAM=FLNM // '.DAT'
  129.     IRSIZE=ICHANS*N*2
  130.     IDELTMS=NINT(DELT*1.0E6)
  131.  
  132. *    Design the bandpass filter.
  133.  
  134.     CALL BPDES (F1,F2,DELT,NS,A,B,C,D,E,GRAF)
  135.     CALL BPDES (F12,F22,DELT,NS,A2,B2,C2,D2,E2,GRAF2)
  136.  
  137. *    Initialize past values in filter to zero.
  138.  
  139.     X=0.0
  140.     X2=0.0
  141.  
  142. *    Write the data in the form of binary numbers to a data
  143. *    file that may be read by SPAD.
  144.  
  145.     OPEN (NUMO,FILE=FILNAM,STATUS='UNKNOWN',ACCESS='SEQUENTIAL',
  146.      +        FORM='BINARY')
  147.     WRITE (NUMO) ICHANS,IRSIZE,NUMREC,IDELTMS
  148.     WRITE (NUMO) (GAIN(I),I=0,7)
  149.     
  150. *    Put message on screen.
  151.  
  152.     WRITE (*,'(/////////////////////16X,
  153.      +           ''C O L O R E D   N O I S E   U T I L I T Y'')')
  154.     WRITE (*,'(/16X,
  155.      +           ''         T W O  C H A N N E L S          '')')
  156.     WRITE (*,'(/25X,''Creating '',A,'' now.''/)') FILNAM 
  157.  
  158. *    Create NUMREC records of the sequence.
  159.  
  160.     DO IB=1,NUMREC
  161.  
  162.       IF (IB .EQ. 1) ITOSS=1000
  163.       IF (IB .NE. 1) ITOSS=0
  164.  
  165. *      Generate the sequence NDATA(J).
  166.  
  167.       DO J=1,N+ITOSS
  168.  
  169.         X(1,5) =SQRT(6.0/DELT)*(RAN1(IDUM)-0.5)
  170.         X2(1,5)=SQRT(6.0/DELT)*(RAN1(IDUM)-0.5)
  171.  
  172. *        Go thru NS sections of filter.
  173.  
  174.         DO I=1,NS
  175.           TEMP=A(I)*(X(I,5)-2.0*X(I,3)+X(I,1))-B(I)*X(I+1,4)
  176.           X(I+1,5)=TEMP-C(I)*X(I+1,3)-D(I)*X(I+1,2)-E(I)*X(I+1,1)
  177.           TEMP2=A2(I)*(X2(I,5)-2.0*X2(I,3)+X2(I,1))-B2(I)*X2(I+1,4)
  178.           X2(I+1,5)=TEMP2-C2(I)*X2(I+1,3)-D2(I)*X2(I+1,2)
  179.      +                  -E2(I)*X2(I+1,1)
  180.         ENDDO
  181.  
  182. *        Push down past values in filter.
  183.  
  184.         DO I=1,NS+1
  185.           DO K=1,4
  186.             X (I,K)=X (I,K+1)
  187.             X2(I,K)=X2(I,K+1)
  188.           ENDDO
  189.         ENDDO
  190.  
  191. *        Throw away ITOSS values to eliminate startup transient.
  192.  
  193.         IF (IB .EQ. 1) THEN
  194.           IF (J .LE. ITOSS) CYCLE
  195.           J2JM1=2*(J-ITOSS)-1
  196.           J2J  =2*(J-ITOSS)
  197.           NDATA(J2JM1)=INT2(X (NS+1,5)/VOFST)
  198.           NDATA(J2J  )=INT2(X2(NS+1,5)/VOFST)
  199.           CYCLE
  200.         ENDIF
  201.  
  202. *        Set NDATA(J) equal to the output of the filter and continue.
  203.  
  204.         J2JM1=2*J-1
  205.         J2J  =2*J
  206.         NDATA(J2JM1)=INT2(X (NS+1,5)/VOFST)
  207.         NDATA(J2J  )=INT2(X2(NS+1,5)/VOFST)
  208.  
  209.       ENDDO
  210.  
  211. *      Display record number message.
  212.  
  213.       IF (IB .EQ. 1) THEN
  214.         WRITE (*,50) IB
  215. 50        FORMAT (25X,'Writing Record ',I4.4)
  216.       ELSE
  217.         WRITE (*,60) IB
  218. 60        FORMAT ('+',24X,'Writing Record ',I4.4)
  219.       ENDIF
  220.  
  221. *      Save output to a file.
  222.  
  223.       WRITE (NUMO) (NDATA(J),J=1,2*N)
  224.  
  225.     ENDDO
  226.  
  227.     CLOSE (NUMO,STATUS='KEEP')
  228.  
  229.     WRITE (*,'( )')
  230.     STOP '                    Program terminated successfully.'
  231.     END
  232.  
  233. *    Bandpass filter routine.   (BPDES.FOR)
  234.  
  235. *    Bandpass Butterworth Digital Filter Design Subroutine
  236.  
  237. *    Inputs are passband (-3 dB) frequencies F1 and F2 in Hz,
  238. *            Sampling interval T in seconds, and
  239. *            number NS of filter sections.
  240.  
  241. *    Outputs are NS sets of filter coefficients, i.e.,
  242. *            A(K) thru E(K) for K = 1 thru NS, and
  243. *            20 pairs of frequency and power gain, i.e.,
  244. *            GRAF(1,K) thru GRAF(2,K) FOR K = 1 thru 20.
  245.  
  246. *    Note that A(K) thru E(K) as well as GRAF(2,20) must be
  247. *            dimensioned in the calling program.
  248.  
  249. *    The digital filter has NS sections in cascade.  The Kth
  250. *            section has the transfer function:
  251.  
  252. *                A(K)*(Z**4-2*Z**2+1)
  253. *        H(Z) =  ------------------------------------
  254. *            Z**4+B(K)*Z**3+C(K)*Z**2+D(K)*Z+E(K)
  255.  
  256. *    Thus, if F(M) and G(M) are the input and output of the 
  257. *            Kth section at time M*T, then
  258.  
  259. *        G(M)=A(K)*(F(M)-2*F(M-2)+F(M-4))-B(K)*G(M-1)
  260. *            -C(K)*G(M-2)-D(K)*G(M-3)-E(K)*G(M-4)
  261.  
  262. *    Do not use this routine for lowpass design.  Use LPDES.
  263. *    Do not use this routine for highpass design.  Use HPDES.
  264.  
  265.     SUBROUTINE BPDES (F1,F2,T,NS,A,B,C,D,E,GRAF)
  266.  
  267.     IMPLICIT     REAL*4 (A-H,O-Z), INTEGER*2 (I-N)
  268.     REAL*4         A(*),B(*),C(*),D(*),E(*),GRAF(2,20)
  269.  
  270. *    PI=2.0*ASIN(1.0)
  271.     PI=3.1415926536
  272.  
  273.     W1=SIN(F1*PI*T)/COS(F1*PI*T)
  274.     W2=SIN(F2*PI*T)/COS(F2*PI*T)
  275.     WC=W2-W1
  276.     Q=WC*WC+2.0*W1*W2
  277.     S=W1*W1*W2*W2
  278.  
  279.     DO 150 K=1,NS
  280.       CS=COS(FLOAT(2*(K+NS)-1)*PI/FLOAT(4*NS))
  281.       P=-2.0*WC*CS
  282.       R=P*W1*W2
  283.       X=1.0+P+Q+R+S
  284.       A(K)=WC*WC/X
  285.       B(K)=(-4.0-2.0*P+2.0*R+4.0*S)/X
  286.       C(K)=(6.0-2.0*Q+6.0*S)/X
  287.       D(K)=(-4.0+2.0*P-2.0*R+4.0*S)/X
  288.       E(K)=(1.0-P+Q-R+S)/X
  289. 150    CONTINUE
  290.  
  291.     DO 170 J=1,2
  292.       DO 160 I=1,10
  293.         K=I*(2-J)+(21-I)*(J-1)
  294.         GRAF(2,K)=0.01+0.98*FLOAT(I-1)/9.0
  295.         X=(1.0/GRAF(2,K)-1.0)**(1.0/FLOAT(4*NS))
  296.         SQ=SQRT(WC*WC*X*X+4.0*W1*W2)
  297.         GRAF(1,K)=ABS(ATAN(0.5*(WC*X+FLOAT(2*J-3)*SQ)))/(PI*T)
  298. 160      CONTINUE
  299. 170    CONTINUE
  300.  
  301.     RETURN
  302.     END
  303.  
  304.     FUNCTION RAN1(IDUM)
  305.  
  306. *    Returns a uniform random deviate between 0.0 and 1.0.  Set
  307. *    IDUM to any negative value to initialize or reinitialize the
  308. *    sequence. Taken from Numerical Recipes, p.196-197.
  309.  
  310.     INTEGER*2    IDUM
  311.     REAL*4        R(97)
  312.     PARAMETER (M1=259200,IA1=7141,IC1=54773,RM1=3.8580247E-6)
  313.     PARAMETER (M2=134456,IA2=8121,IC2=28411,RM2=7.4373773E-6)
  314.     PARAMETER (M3=243000,IA3=4561,IC3=51349)
  315.     DATA        IFF /0/
  316.  
  317.     IF (IDUM.LT.0.OR.IFF.EQ.0) THEN
  318.       IFF=1
  319.       IX1=MOD(IC1-IDUM,M1)
  320.       IX1=MOD(IA1*IX1+IC1,M1)
  321.       IX2=MOD(IX1,M2)
  322.       IX1=MOD(IA1*IX1+IC1,M1)
  323.       IX3=MOD(IX1,M3)
  324.       DO J=1,97
  325.         IX1=MOD(IA1*IX1+IC1,M1)
  326.         IX2=MOD(IA2*IX2+IC2,M2)
  327.         R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1
  328.       ENDDO
  329.       IDUM=1
  330.     ENDIF
  331.  
  332.     IX1=MOD(IA1*IX1+IC1,M1)
  333.     IX2=MOD(IA2*IX2+IC2,M2)
  334.     IX3=MOD(IA3*IX3+IC3,M3)
  335.     J=1+(97*IX3)/M3
  336.     IF(J.GT.97.OR.J.LT.1)PAUSE
  337.     RAN1=R(J)
  338.     R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1
  339.  
  340.     RETURN
  341.     END
  342.